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constraints  are  satisfied  to  first  order.  During  this 
iteration,  the  first  variation  of  the  functional  is  minimized, 
subject  to  the  linearized  constraints  and  to  a quadratic 
constraint  imposed  on  the  variations  of  the  control,  the 
parameter,  and  the  missing  components  of  the  initial  state. 

The  restoration  phase  involves  one  or  more  iterations  and  is 
designed  to  force  constraint  satisfaction  to  a predetermined 
accuracy,  while  the  norm  squared  of  the  variations  of  the 
control,  the  parameter,  and  the  missing  components  of  the 
initial  state  is  minimized.  ^ — — 

The  algorithm  has  two  principal  properties : (i)  it 
produces  a sequence  of  feasible  suboptimal  solutions;  the 
functions  obtained  at  the  end  of  each  cycle  satisfy  the 
constraints  to  a predetermined  accuracy;  and  (ii)  the  value 
of  the  functional  at  the  end  of  any  complete  conjugate 
gradient-restoration  cycle  is  smaller  than  the  value  of  the 
same  functional  at  the  beginning  of  that  cycle. 

The  sequential  conjugate  gradient-restoration  algorithm 
presented  here  differs  from  previous  work,  in  that  it  is  not 
required  that  the  state  vector  be  given  at  the  initial  point. 
Instead,  the  initial  conditions  can  be  absolutely  general. 
Since  the  present  algorithm  is  capable  of  handling  general 
final  conditions,  it  is  suitable  for  the  solution  of  optimal 
control  problems  with  general  boundary  conditions.  Its 
importance  lies  in  the  fact  that  many  optimal  control  problems 
involve  initial  conditions  of  the  type  considered  here. 

Nine  numerical  examples  are  presented  in  order  to  illus- 
trate the  performance  of  the  algorithm.  The  numerical 
results  show  the  feasibility  as  well  as  the  convergence 
characteristics  of  the  present  algorithm. 
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Sequential  Conjugate  Gradient-Restoration  Algorithm 
for  Optimal  Control  Problems 
with  General  Boundary  Conditions^ 


A.K.  and  A.  MIELE^ 


Abstract.  This  paper  considers  the  numerical  solution  of 
the  problem  of  minimizing  a functional  I subject  to  differen- 
tial constraints  and  general  boundary  conditions.  It  consists 
of  finding  the  state  x{t),  the  control  u(t),  and  the  parameter 
IT  so  that  the  functional  I is  minimized  while  the  constraints 
and  the  boundary  conditions  are  satisfied  to  a predetermined 
accuracy. 
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The  approach  taken  is  a sequence  of  two-phase  cycles, 
composed  of  a conjugate  gradient  phase  and  a restoration 
phase.  The  conjugate  gradient  phase  involves  one  iteration 
and  is  designed  to  decrease  the  value  of  the  functional,  while 
the  constraints  are  satisfied  to  first  order.  During  this 
iteration,  the  first  variation  of  the  functional  is  minimized, 
subject  to  the  linearized  constraints.  The  minimization  is 
performed  over  the  class  of  variations  of  the  control,  the 
parameter,  and  the  missing  components  of  the  initial  state 
which  are  equidistant  from  some  constant  multiple  of  the  cor- 
responding variations  of  the  previous  conjugate  gradient  phase. 
The  sequence  of  conjugate  gradient  phases  generated  by  the 
algorithm  is  such  that,  for  the  special  case  of  a quadratic 
functional  subject  to  linear  constraints,  various  orthogona- 
lity and  conjugacy  conditions  hold.  The  restoration  phase  in- 
volves one  or  more  iterations  and  is  designed  to  force  cons- 
traint satisfaction  to  a predetermined  accuracy,  while  the  norm 
squared  of  the  variations  of  the  control,  the  parameter,  and 
the  missing  components  of  the  initial  state  is  minimized. 

The  principal  property  of  the  algorithm  is  that  it  pro- 
duces a sequence  of  feasible  suboptimal  solutions;  the  fun- 
ctions obtained  at  the  end  of  each  cycle  satisfy  the  cons- 
traints to  a predetermined  accuracy.  Therefore,  the  values 
of  the  functional  I corresponding  to  any  two  elements  of  the 
sequence  are  comparable. 
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The  stepsize  of  the  conjugate  gradient  phase  is  deter- 
mined by  a one-dimensional  search  on  the  augmented  functional 
J,  while  the  stepsize  of  the  restoration  phase  is  obtained  by 
a one-dimensional  search  on  the  constraint  error  P.  The  con- 
jugate gradient  stepsize  and  the  restoration  stepsize  are  chosen 
so  that  the  restoration  phase  preserves  the  descent  property  of  the 
conjugate  gradient  phase.  Therefore,  the  value  of  the  functional  I 
at  the  end  of  any  complete  conjugate  gradient-restoration 
cycle  is  smaller  than  the  value  of  the  same  functional  at  the 
beginning  of  that  cycle.  Of  course,  restarting  the  algorithm 
might  be  occasionally  necessary. 

The  sequential  conjugate  gradient-restoration  algorithm 
presented  here  differs  from  that  of  Refs.  3 and  4,  in  that  it 
is  not  required  that  the  state  vector  be  given  at  the  initial 
point.  Instead,  the  initial  conditions  can  be  absolutely 
general.  In  analogy  with  Refs.  3 and  4,  the  present  algorithm 
is  capable  of  handling  general  final  conditions;  therefore, 
it  is  suitable  for  the  solution  of  optimal  control  problems 
with  general  boundary  conditions.  Its  importance  lies  in  the 
fact  that  many  optimal  control  problems  involve  initial  con- 
ditions of  the  type  considered  here. 

Nine  numerical  examples  are  presented  in  order  to  illus- 
trate the  performance  of  the  algorithm.  The  numerical  results 
show  the  feasibility  as  well  as  the  convergence  characteris- 
tics of  the  present  algorithm. 
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Key  Words.  Optimal  control,  numerical  methods,  computing 
methods,  gradient  methods,  gradient-restoration  algorithms, 
sequential  gradient-restoration  algorithms,  conjugate  gradient- 
restoration  algorithms,  sequential  conjugate  gradient- 
restoration  algorithms,  general  boundary  conditions. 
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1 . Introduction 

Over  the  past  several  years,  a successful  family  of  al- 
gorithms for  the  solution  of  optimal  control  problems  involv- 
ing differential  constraints  and  terminal  constraints  has  been 
developed  at  Rice  University  by  Miele  and  his  associates  (see 
Refs.  1-4) . They  are  known  as  sequential  gradient-restoration 
algorithms  and  have  been  designed  for  the  solution  of  differ- 
ent classes  of  optimal  control  problems.  Some  of  these  algo- 
rithms are  of  the  ordinary-gradient  type  (Refs.  1-2),  while 
the  rest  are  of  the  conjugate-gradient  type  (Refs.  2-4).  All 
of  them  have  shown  to  be  robust  and  reliable;  however,  they 
all  require  the  state  vector  to  be  given  at  the  initial  point. 
Owing  to  the  fact  that  optimal  control  problems  exist,  which 
require  satisfaction  of  more  general  boundary  conditions,  the 
task  of  extending  the  aforementioned  family  of  algorithms  must 
be  undertaken,  with  these  ideas  in  mind;  (i)  to  retain  the 
robustness,  reliability,  and  convergence  characteristics  of 
the  algorithms  discussed  in  Refs.  1-4;  (ii)  to  be  able  to 
handle  all  of  the  optimal  control  problems  treated  in  Refs. 
1-4;  and  (iii)  to  have  the  additional  capability  of  handling 
optimal  control  problems  with  general  boundary  conditions. 

In  Ref.  5,  an  algorithm  of  the  ordinary  gradient  type 
was  developed,  extending  the  work  of  Ref.  1 to  problems  with 
general  boundary  conditions.  In  this  report,  an  algorithm  of 
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the  conjugate  gradient  type  is  developed,  extending  the  work 
of  Refs.  3-4  to  problems  withgeneral  boundary  conditions. 

Specifically,  the  following  optimal  control  problem  is 
considered:  Minimize  the  functional  I,  which  depends  on  the 

n-vector  state  x(t),  the  m-vector  control  u(t),  and  the 
p-vector  parameter  ir.  The  state  and  the  parameter  are  re- 
quired to  satisfy  r scalar  relations  at  the  initial  point  and 
q scalar  relations  at  the  final  point.  Along  the  interval  of 
integration,  the  state,  the  control,  and  the  parameter  are 
required  to  satisfy  n scalar  differential  equations. 

The  approach  taken  is  a sequence  of  two-phase  cycles, 
composed  of  a conjugate  gradient  phase  and  a restoration  phase. 
The  conjugate  gradient  phase  involves  one  iteration  and  is  de- 
signed to  decrease  the  value  of  the  functional,  while  the  con- 
straints are  satisfied  to  first  order.  During  this  iteration, 
the  first  variation  of  the  functional  is  minimized,  subject  to 
the  linearized  constraints.  The  minimization  is  performed 
over  the  class  of  variations  of  the  control,  the  parameter, 
and  the  missing  components  of  the  initial  state  which  are 
equidistant  from  some  constant  multiple  of  the  corresponding 
variations  of  the  previous  conjugate  gradient  phase.  The 
sequence  of  conjugate  gradient  phases  generated  by  the  algo- 
rithm is  such  that,  for  the  special  case  of  a quadratic 
functional  subject  to  linear  constraints,  various  orthogona- 
lity and  conjugacy  conditions  hold.  The  restoration  phase 
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involves  one  or  more  iterations  and  is  designed  to  force  cons- 
traint satisfaction  to  a predetermined  accuracy, while  the  norm 
squared  of  the  variations  of  the  control,  the  parameter,  and 
the  missing  components  of  the  initial  state  is  minimized. 

The  principal  property  of  the  algorithm  is  that  it  pro- 
duces a sequence  of  feasible  suboptimal  solutions;  the  fun- 
ctions obtained  at  the  end  of  each  cycle  satisfy  the  cons- 
traints to  a predetermined  accuracy.  Therefore,  the  values  of 
the  functional  I corresponding  to  any  two  elements  of  the  se- 
quence are  comparable. 

The  stepsize  of  the  conjugate  gradient  phase  is  deter- 
mined by  a one-dimensional  search  on  the  augmented  functional 
J,  while  the  stepsize  of  the  restoration  phase  is  obtained  by 
a one-dimensional  search  on  the  constraint  error  P.  The  con- 
jugate gradient  stepsize  and  the  restoration  stepsize  are  chosen 
so  that  the  restoration  phase  preserves  the  descent  property  of  the 
conjugate  gradient  phase.  Therefore,  the  value  of  the  functional  I 
at  the  end  of  any  complete  conjugate  gradient-restoration 
cycle  is  smaller  than  the  value  of  the  sane  functional  at  the 
beginning  of  that  cycle.  Of  course,  restarting  the  algorithm 
might  be  occasionally  necessary. 

A time  normalization  is  used  in  order  to  simplify  the 
numerical  computations.  Specifically,  the  actual  time  0 is 
replaced  by  the  normalized  time  t = 0/t , which  is  defined  in 
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such  a way  that  the  initial  time  is  t = 0 and  the  final  time 
is  t=  1.  The  actual  final  time  T,  if  it  is  free,  is  regarded 
as  a component  of  the  vector  parameter  tt  to  be  optimized.  In 
this  way,  an  optimal  control  problem  with  variable  final  time 
is  converted  into  an  optimal  control  problem  with  fixed  final 
time. 
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2.  Statement  of  the  Problem 

2.1.  Notation.  Let  t denote  the  independent  variable, 
and  let  x(t)  , u(t),  it  denote  the  dependent  variables.  The 
time  t is  a scalar,  the  state  x(t)  is  an  n-vector,  the  con- 
trol u(t)  is  an  m-vector,  and  the  parameter  tt  is  a p-vector,^ 

The  state  x(t)  is  partitioned  into  vectors  y(t)  and  z(t), 
defined  as  follows;  y(t)  is  an  a-vector  including  those  com- 
ponents of  the  state  that  are  prescribed  at  the  initial  point, 
and  z(t)  is  a b-vector  including  those  components  of  the 
state  that  are  not  prescribed  at  the  initial  point.  Clearly, 
a + b = n. 

2.2.  Optimization  Problem.  With  the  above  notations, 
the  optimization  problem  can  be  stated  as  follows.  Minimize 
the  functional 


1=  I f (x,u,TT,  t)dt  + [h(z,Tr)]^  + [g(x,Tr)],  , 


with  respect  to  the  state  x(t),  the  control  u(t),  and  the 
parameter  tt  which  satisfy  the  differential  constraints 


X - (f)  (X,U,  TT,  t)  - 0 , 


and  the  boundary  conditions 


All  vectors  are  column  vectors. 


0 < t < 1, 
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y (0)  = given. 

(3) 

[w(z,Tr)  }q  = 0 , 

(4) 

I'/'  (x,7r)  ] = 0 . 

(5) 

In  the  above  equations,  the  quantities 

I,f,h,g  are  scalar,  the 

function  ((>  is  an  n-vector,  the  function  u is  a c-vector,  and 
the  function  is  a q-vector.  The  number  of  initial  condi- 
tions r = a + c and  the  number  of  final  conditions  q must  satis- 
fy the  following  relation: 

r + q < n^  + n^^  + p*  < 2n  + p . (6) 

In  Ineq.  (6) , the  symbol  p*  < P denotes  the  number  of  compo- 
nents of  the  parameter  ir  present  in  the  boundary  conditions; 
the  symbol  ng  < n denotes  the  number  of  state  variables  pres- 
ent in  the  initial  conditions;  and  the  symbol  n^^  < n denotes 
the  number  of  state  variables  present  in  the  final  conditions. 

2.3.  First-Order  Conditions.  From  calculus  of  varia- 
tions, it  can  be  seen  that  the  previous  problem  is  one  of  the 
Bolza  type,  and  it  can  be  recast  as  that  of  minimizing  the 
augmented  functional 


J=I  + L,  (7-1) 

subject  to  (2)-(5).  In  compact  notation,  the  functional  I 
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can  be  rewritten  as 


= \ fdt+  (h)Q+  (g)^  , 
Jn 


(7-2) 


and  the  functional  L is  defined  as 


L = \ x'^(x  - ((i)dt  + (a'^o))  Q + , 


(7-3) 


' h? 
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where  the  n-vector  X (t)  is  a variable  Lagrange  multiplier, 
the  c-vector  a is  a constant  Lagrange  multiplier,  and  the 
q-vector  y is  a constant  Lagrange  multiplier.  After  integra- 
ting by  parts  the  term  x'^x,  the  functional  (7-3)  can  be 
rewritten  as 


\ (“X'^x  - x'^((i)dt  + (-\'^x+ 0*^(0)  Q + (X'^x  + y’^t/;)  . (7-4) 


The  functions  x(t),  u(t),  it  and  the  multipliers  X (t)  , a, 
y solving  the  previous  problem  must  satisfy  the  feasibility 
equations  (2) -(5)  and  the  following  optimality  conditions; 


X - fx  + <^x^  = 0, 


f - ())  X = 0, 
u u 


0 < t < 1, 


0 < t < 1, 


(8) 

(9) 


(f^-  <(>^x)dt+  (h^  + w^a)(j+  (g^  + 'P^y)i  = 0,  (10) 
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{-C+ + a)^a)Q  = 0,  (11) 

(X  + gx  + ’l'x^)l  = 0 . (12) 

2.4.  Remark . Just  as  the  state  vector  x(t)  is  parti- 

tioned into  an  a-vector  y(t)  and  a b-vector  z(t),  the  multi- 
plier vector  \(t)  is  partitioned  into  a a-vector  n(t)  and  a 
b-vector  ^(t) , having  the  following  meaning:  n (t)  is  associa- 
ted with  y(t),  and  ^(t)  is  associated  with  2(t).  With 
reference  to  Eq.  (11) , C(0)  denotes  the  portion  of  the  initial 

Lagrange  multiplier  X(0)  which  is  associated  with  z(0),  the 
portion  of  the  initial  state  vector  x(0)  which  is  not  pres- 
cribed. 

2.5.  Approximate  Methods.  Since  in  general  the  differ- 
ential system  (2)- (5)  and  (8)- (12)  is  nonlinear,  approximate 
methods  are  employed  to  find  a solution  iteratively.  In  this 
connection,  let  the  norm  squared  of  a vector  v be  defined  by 

N(v)  = v^v  . (13) 

Then,  the  constraint  error  P can  be  written  as^ 

P=1  Ntx- (}))dt  + N(a))Q  + N(il;)j^  , (14) 

Jq 

C 

In  Eq.  (14),  it  is  assumed  that  the  initial  condition  (3)  is 
satisfied. 
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and  the  error  in  the  optimality  conditions  Q is  given  by 


Q = 


\ N(X  - 

Jo 


f + (t,  x)dt  + 
X X 


I N(f  - (|) 

I u 

Jn 


X)dt 


r ri 


+ N 


(f^-4.^X)dt+  (h^  + (.^a)o+(g^+,|i^M)^ 


+ N (-C  + h^  + 0 N(X  + g^  + • 


(15) 


For  the  exact  optimal  solution,  one  must  have 


P=0,  Q=0.  (16) 

For  an  approximation  to  the  optimal  solution,  the  following 
relations  are  to  be  satisfied; 
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3 . Description  of  the  Algorithm 

The  technique  employed  is  characterized  by  a sequence  of 
two-phase  cycles,  composed  of  a conjugate  gradient  phase  and 
a restoration  phase.  These  phases  are  described  below. 

The  conjugate  gradient  phase  is  started  only  when  the 
constraint  error  P satisfies  Ineq.  (17-1) . It  involves  a 
single  iteration,  which  is  designed  to  decrease  the  value  of 
the  functional  I or  the  augmented  functional  J,  while  the 
constraints  are  satisfied  to  first  order.  During  this  itera- 
tion, the  first  variation  of  the  functional  I is  minimized, 
subject  to  the  linearized  constraints.  The  minimization  is 
performed  over  the  class  of  variations  of  the  control  u(t)  and 
the  parcuneters  it  and  z(0)  which  are  equidistant  from  some  con- 
stant multiple  of  the  corresponding  variations  of  the  previous 
conjugate  gradient  phase. ^ 

The  restoration  phase  is  started  only  when  the  constraint 
error  P violates  Ineq.  (17-1) . The  restoration  phase  involves 
one  or  more  iterations.  In  each  restorative  iteration,  the 
objective  is  to  reduce  the  constraint  error  P,  while  the 


The  sequence  of  conjugate  gradient  phases  generated  by  the 
algorithm  is  such  that,  for  the  special  case  of  a quadratic 
functional  subject  to  linear  constraints,  various  orthogonal- 
ity and  conjugacy  conditions  hold  (see  Section  6 ) . 
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constraints  are  satisfied  to  first  order  and  the  norm  squared 
of  the  variations  of  the  control  u(t)  and  the  parameters  tt 
and  z(0)  is  minimized.  The  restoration  phase  is  terminated 
whenever  Ineq.  (17-1)  is  satisfied. 

A complete  conjugate  gradient-restoration  cycle  is  de- 
signed so  that  the  value  of  the  functional  I decreases,  while 
the  constraints  are  satisfied  to  the  accuracy  (17-1)  both  at 
the  beginning  and  at  the  end  of  the  cycle.  Finally,  the  algo- 
rithm as  a whole  is  terminated  whenever  Ineqs.  (17)  are  satis- 
fied simultaneously. 

3.1.  Remark.  During  each  conjugate  gradient  iteration 
or  restorative  iteration,  use  of  nonlinear  equations  must  be 
avoided.  Therefore,  the  exact  feasibility  equations  (2) -(5) 
are  replaced  by  their  corresponding  linearized  approximations. 
These  linearized  approximations  do  not  include  forcing  terms 
in  the  conjugate  gradient  phase,  while  they  do  include  forcing 
terms  in  the  restoration  phase. 

3.2.  Notation.  For  any  iteration  of  the  conjugate 
gradient  phase  or  the  restoration  phase,  the  following  termi- 
nology is  adopted:  x(t),  u(t),  it  denote  the  nominal  functions; 
x(t),  u(t),  TT  denote  the  varied  functions;  and  Ax(t),  Au(t),  Air 
denote  the  displacements  leading  from  the  nominal  functions  to 
the  varied  functions.  These  quantities  satisfy  the  relations 

X(t)=X(t)+Ax(t),  U(t)=u(t)+Au(t),  Tf=TT  + AlT  . (18) 


f 
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Because  the  vector  x is  partitioned  into  y and  z,  Eq.  (18-1) 
implies  that 

y (t)  = y (t)  + Ay (t) , z(t)  = z(t)  + Az(t) . (19) 

Let  a be  a positive  number  representing  the  stepsize 
(either  the  gradient  stepsize  or  the  restoration  stepsize) . 
Then,  we  define  the  displacements  per  unit  stepsize  as  fol- 
lows : 

A(t)  = Ax(t)/a,  B(t)  = Au(t)/a,  C=ATr/a.  (20) 

The  vector  A is  partitioned  into  vectors  D and  E associated 
with  y and  z,  respectively.  Therefore,  Eq.  (20-1)  implies 
that 

D(t)  = Ay (t)/a,  E(t)  = Az (t)/a.  (21) 

Upon  combining  (18)- (19)  with  (20) -(21),  we  see  that 

x(t)  = x(t)  + aA(t)  , u (t)  = u (t)  + aB  (t)  , rf  = tt  + aC,  (22) 

y (t)  = y (t)  + aD(t)  , z (t)  = z (t)  + aE (t)  . (23) 

3*3.  Desired  Properties.  The  functions  Ax(t),  Au(t), 

Att  must  be  determined  so  as  to  produce  some  desirable  effect 
at  every  iteration,  namely,  the  decrease  of  the  functionals  I, 
and/or  J,  and/or  P.  Thus,  the  following  descent  properties 
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I < I,  and/or  J < J,  and/or  P < F, 


where  I,  J,  P are  associated  with  the  nominal  functions  and  I, 
J,  P are  associated  with  the  varied  functions.  In  turn,  the 
functions  A(t) , B(t),  C are  chosen  so  that 


61 < 0,  and/or  6J  < 0,  and/or  6P  < 0, 


where  the  symbol  6(...)  denotes  the  first  variation.  Then,  by 
choosing  the  stepsize  a sufficiently  small,  the  satisfaction 
of  relations  (24)  is  guaranteed.  Ineqs.  (24-1) , (24-2)  and 
(25-1),  (25-2)  characterize  the  conjugate  gradient  phase,  while 
Ineqs.  (24-3)  and  (25-3)  characterize  the  restoration  phase. 

3.4.  First  Variations.  Next,  we  give  the  expressions 
for  the  first  variations  of  the  functionals  I,  J,  P;  after 
simple  manipulations,  omitted  for  the  sake  of  brevity,  they 
take  the  form^'® 


6l/a=  \ (f^A+f^B+f^C)dt+  (h^E  + hJc)o+  (gJ^A  + gJJc)^, 


^Implicit  in  Eqs.  (26)-(28)  is  the  assumption  D(0)  =0. 

®The  first  variation  of  the  augmented  functional  J is  computed 
by  varying  the  functions  x(t),  u(t),  tt,  while  holding  the 
multipliers  A(t),o,u  unchanged. 
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and 


6J/a=  \ (-A  + f^- ((.^X)  ^Adt  + 


r 

1 (f  -4> 

1 u 

Jo 


X ) "^Bdt 


r 

\ (f  - (j) 

1 ' IT  ^TT 

Jo 


IT 


X)dt+  (h^  + u)^a)Q+  + 


+ [^  (-C +h^  + a)^a)'^E]^+  [(^+ 9^  + ' 


(27 


and 


6P/2a=  \ (X- 4))'^(A- (t)^A- ({)Jb- (|)JJc)dt 


+ [^  0)^(0) ^E  + (i) Jo  + [^i^^(ipjA+ti-Jc)j  . 


T T 

I'  + 


(28 


3.5.  Remark.  For  the  purposes  of  this  report,  Eqs, 

(26) -(28)  must  be  completed  by  one  of  the  following  relations 


K/a^  = \ B'^Bdt  + c'^C  + (e'^E)  q , 


(29; 


or 


( (B-yB)'*’(B-^ 


K/a‘=  \ (B-YB)^(B-YB)dt+(C-YC)"(C-YC)  + [(E-YE)'^(E-YE)l  , (30 


! f i’: 

■ ■; 


i ! . 
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which  constitute  a measure  of  the  overall  change  of  the  con- 
trol u(t)  and  the  parameters  ir  and  z(0).  Equation  (29)  is 
employed  in  the  restoration  phase,  and  Eq.  (30)  is  employed  in 
the  conjugate  gradient  phase.  In  Eq.  (30) , the  functions 
A(t),  B(t),  C pertain  to  the  present  conjugate  gradient  phase, 
and  the  functions  A(t) , B(t),  C pertain  to  the  previous  con- 
jugate gradient  phase.  The  symbol  y denotes  a scalar,  non- 
negative quantity,  called  the  directional  coefficient.  Its 
specification  is  discussed  in  Section  6. 
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4 . Restoration  Phase 

4.1.  Linearized  Equations.  Let  x(t)  , u{t),  tt  denote 
nominal  functions  satisfying  (3),  but  not  necessarily  (2),  (4), 
(5) . To  first  order,  the  perturbations  per  unit  stepsize 
A(t),  B(t),  C must  satisfy  the  linearized  constraint  equations 


A - (f)3A-  4)3c+  (x-  4))  = 0, 

A U TT 

0 < t < 1, 

(31) 

D(0)  = 0, 

(32) 

(dJ^E  + lii^C  + co)  = 0 / 

Z TT  0 

(33) 

(4»^a  + 4;JJc  + 4^)3^  = 0 . 

(34) 

4.2.  Descent  Property.  The  linearized  equations 
(31) -(34)  admit  an  infinite  number  of  solutions,  each  of  which 
is  characterized  by  a descent  property  in  the  constraint  error 
P.  This  descent  property  can  be  seen  by  combining  (28)  with 
(31)-(34);  the  first  variation  of  P becomes 


5P  = -2aP.  (35) 

Since  P>  0,  Eq.  (35)  shows  that  5P < 0;  hence,  for  a suffi- 
ciently small,  a decrease  in  the  constraint  error  P is  guar- 
anteed . 

4.3.  Special  Variations.  Among  the  infinite  number  of 
solutions  of  Eqs.  (31)-(34),  the  one  that  minimizes  the 
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functional  (29)  is  selected.  Thus,  we  seek  the  minimum  of  the 
quadratic  functional  (29),  with  respect  to  the  functions  A(t) , 
B(t),  C which  satisfy  the  linearized  constraints  (31)-(34). 

By  applying  standard  techniques  of  optimal  control  theory 
or  calculus  of  variations,  the  following  optimality  conditions 
are  obtained: 


B=(|)^X,  0<t<l,  (36) 

C=\  <^^Xdt-  (a)^a)Q-  (il^^y)^^,  (37) 

Jo 

E(0)  = (c  - (38) 

i + (t.j^X  = 0,  0<t<l,  (39) 

(^  + = 0.  (40) 


Summarizing,  we  seek  functions  A(t) , B(t),  C and  multipliers 
X (t) , a,  u which  satisfy  the  linearized  constraints  (31) -(34) 
and  the  optimality  conditions  (36)-(40). 

4.4.  Linear,  Two-Point  Boundary-Value  Problem.  The 
technique  used  to  solve  the  LTP-BVP  (31) -(34)  and  (36) -(40) 
is  a backward- forward  integration  scheme  in  combination  with 
the  method  of  particular  solutions  (Refs. 6-8) . The  technique 
requires  the  execution  of  q + 1 independent  sweeps  of  the 
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differential  system  (31) -(34)  and  (36) -(40),  each  character- 
ized by  a different  value  of  the  multiplier  y . 

The  generic  sweep  is  started  by  assigning  particular 
values  to  the  components  of  y;  then,  the  multiplier  A(l)  is 
obtained  from  (40) . Next,  Eq.  (39)  is  integrated  backward  to 
obtain  the  function  X (t) , and  Eq.  (36)  is  employed  to  obtain 
B(t).  With  X(0)  known,  C(0)  is  known.  Therefore,  Eqs.  (33), 
(37),  (38)  constitute  a system  of  b + c + p linear  relations  in 
which  the  unknowns  are  the  b + c + p components  of  the  vectors 
E(0),  a,  C.  For  this  system  to  have  a unique  solution,  the 

9 

following  disequation  must  hold; 

det  1 +a)'‘aj  I ^ 0.  (41) 

L z z TT  ttJo 

With  E(0)  known  and  because  of  (32),  the  vector  A(0)  is  known. 
Then,  A(t)  is  obtained  by  forward  integration  of  (31)  . In 
this  way,  the  sweep  is  completed;  for  the  arbitrary  value 
assigned  to  y , it  leads  to  the  satisfaction  of  all  of  the 
equations  of  the  system  (31)-(34)  and  (36)-(40),  except  Eq. 
(34)  . 


4 


9 

Disequation  (41)  is  obtained  from  (33) , (37) , (38)  after 

elimination  of  E(0)  and  C.  The  resulting  linear  equation  in 
a admits  a unique  solution  providing  (41)  is  satisfied. 


t 
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In  order  to  satisfy  Eq.  (34)  and  because  the  system 
(31)-(34)  and  (36)-(40)  is  nonhomogeneous , q+1  independent 
sweeps  must  be  executed  employing  q+1  different  multiplier 
vectors  y^,  i = l,...,q  + l.  The  first  q sweeps  are  performed 
by  choosing  the  vectors  y,, ,y  to  be  the  columns  of  the 

X 

identity  matrix  of  order  q.  The  last  sweep  is  executed  by 
choosing  null  vector.  As  a result,  one  gener- 

ates the  functions  and  multipliers 

A^(t),  B^(t),  C^,  A^(t),  0^,  y^,  i=l, ,q+l.  (42) 

Now,  we  introduce  the  q+1  undetermined,  scalar  constants 
and  form  the  linear  combinations 


A(t)  = Ek^A^(t) , 


X(t)  = Ek^A^(t) , 


B(t)  = Ek^B^(t) , 


a = Ek . a . , 
1 1 ' 


C=Ek^C^,  (43) 


y = Ek^y^,  (44) 


where  the  summations  are  taken  over  the  index  i . The  q+1 
coefficients  k^  are  obtained  by  forcing  the  linear  combinations 
(43)  to  satisfy  Eq.  (34)  , together  with  the  normalization  con- 
dition (Ref.  6) 


Ek^  = 1 . 


Once  the  constants  k^  are  known,  the  solution  of  the  LTP-BVP 
(31)-(34)  and  (36)-(40)  is  given  by  (43)-(44). 
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4.5.  Restoration  Stepsize.  With  the  functions  A(t) , 


B(t),  C known,  the  one-parameter  family  of  varied  functions 


(22)  can  be  formed.  For  this  one-parameter  family,  the  con- 


straint error  (14)  becomes  a function  of  the  form 


P = P(a)  . 


Then,  the  stepsize  a must  be  selected  so  that  the  following 


relations  are  satisfied; 


P(a)  < P(0)  , 


T (a)  > 0 . 


Satisfaction  of  Ineq.  (47-1)  is  possible  because  of  the  des- 


cent property  of  the  restoration  phase.  Ineq.  (47-2)  is 


required  for  problems  with  free  final  time. 


In  order  to  achieve  satisfaction  of  (47),  a bisection 


process  is  applied  to  the  restoration  stepsize  a,  starting 


from  the  reference  stepsize  aQ=l.  This  reference  stepsize 


has  the  property  of  yielding  one-step  restoration  for  the  case 


where  the  constraints  (2) -(5)  are  linear. 


4.6.  Iterative  Procedure  for  the  Restoration  Phase.  The 


descent  property  (35)  of  the  restoration  phase  guarantees  sat- 


isfaction of  Ineq.  (47-1)  at  the  end  of  any  iteration,  but 


not  s».  tisfaction  of  Ineq.  (17-1).  Therefore,  the  restoration 


algorithm  must  be  employed  iteratively  until  Ineq.  (17-1)  is 


satisfied.  At  this  point,  the  restoration  phase  is  terminated. 
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5.  Conjugate  Gradient  Phase;  General  Case 

5.1.  Linearized  Equations.  Suppose  that  nominal  func- 
tions x(t),  u(t),  ir  are  available,  which  satisfy  (2) -(5).  To 
first  order,  the  perturbations  per  unit  stepsize  A(t) , B(t), 

C must  satisfy  the  linearized  constraint  equations 


, fp  fp  ^p 

A - (^^A-  4.^C=  0, 


D(0)  = 0, 


0 < t < 1, 


(■o^E  + »Jc)„  = 0, 


5.2.  Special  Variations.  Among  the  infinite  number  of 
solutions  of  Eqs.  (48)-(51),  the  one  that  minimizes  the  func- 
tional (26)  is  selected.  Thus,  we  seek  the  minimum  of  the 
linear  functional  (26),  with  respect  to  the  functions  A(t) , 
B(t),  C which  satisfy  the  linearized  constraints  (48)-(51) 
and  the  quadratic  isoperimetric  constraint  (30) . 

By  applying  standard  techniques  of  optimal  control  theory 
or  calculus  of  variations,  the  following  optimality  conditions 
are  obtained 


In  Eqs.  (52)-(56),  the  multiplier  v associated  with  the 
isoperimetric  constraint  (30)  is  set  at  the  level  2v  = 1 . 
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B = yB-  (f^-  ({.^X)  , 0 < t < 1, 


C = yC-  I (f^- (|)^X)dt+  (h^+a)^0)Q  + + , (53 


E(0)  = YE(0)-(-^+h^+  w^a)Q  , 


^ = 0,  0 < t < 1 , 


(A  + = 0. 


Svunmarizing,  for  a given  value  of  the  directional  coefficient 
y,  we  seek  functions  A(t) , B(t),  C and  multipliers  X (t) , o,  y 
which  satisfy  the  linearized  constraints  (48) -(51)  and  the 
optimality  conditions  (52) -(56). 

5.3.  Isoperimetric  Constraint.  In  the  light  of  (52) -(56) 
the  error  in  the  optimality  conditions (15)  reduces  to 


Q=\  (B-yB)  ^(B-YB)dt  + (C-yC)  ^(C-YC)+  [(E-yE)^  (E-yE)  ] q • (57 


Consequently,  the  following  relation  ties  the  isoperimetric 
constant,  the  stepsize,  and  the  error  in  the  optimality  con- 
ditions : 


K = a Q . 
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Clearly,  to  assign  values  to  the  isoperimetric  constant 
is  the  Scun\e  as  assigning  values  to  the  stepsize.  However, 
there  is  no  clear-cut  way  of  determining  a priori  convenient 
values  for  the  isoperimetric  constant.  Therefore,  the  imple- 
mentation of  the  conjugate  gradient  algorithm  becomes  simpler 
if  one  avoids  evaluating  a in  terms  of  K through  (58)  and 
assigns  values  to  a directly. 

5.4.  Descent  Property.  When  the  variations  defined  by 
(48)- (56)  are  employed,  the  first  variation  of  the  augmented 
functional  (27)  becomes 

6J=-a(Q  + YZ),  (59) 

where  Q is  given  by  (57)  and 

Z=  ( (B- YB)'^Bdt  + (C-YC)'^C+[(E- yE)'^e]q  . (60) 

Jo 

For  the  first  iteration  of  the  conjugate  gradient  phase, 
one  sets 


Y=  0 , 


(61) 


with  the  implication  that 


6J  = -aQ  . 


(62) 
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Since  Q>0,  Eq.  (62)  shows  that  6J  < 0.  Hence,  for  a suffi- 
ciently small,  it  is  guaranteed  that  the  augmented  functional 
J decreases. 

For  subsequent  iterations,  one  sets  Y?^0.  More  specifi- 
cally, the  directional  coefficient  must  be  such  that 

Y>0,  (63) 

and  its  proper  value  is  discussed  in  Section  6.  At  any  rate, 
Eq.  (59)  shows  that  6J  < Q providing 

Q + yZ  > 0 . (64) 

Hence,  for  a sufficiently  small,  it  is  guaranteed  that  the 
augmented  functional  J decreases  as  long  as  Ineq.  (64)  is 
satisfied.  If  Ineq.  (64)  is  violated,  the  descent  property  on 
J no  longer  holds,  and  the  conjugate  gradient  phase  must  be 
restarted  by  resetting  the  directional  coefficient  y at  the 
level  (61). 

5.5.  Linear,  Two-Point  Boundary-Value  Problem.  For  a 
given  value  of  the  directional  coefficient  y,  the  technique 
used  to  solve  the  LTP-BVP  (48) -(56),  associated  with  the  con- 
jugate gradient  phase,  is  analogous  to  that  described  for  the 
restoration  phase  (see  Section  4.4);  hence,  it  is  not  re- 
peated, for  the  sake  of  brevity. 
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5.6.  General  Solution.  Next,  assume  that  two  particular 
values  are  given  to  the  directional  coefficient  y,  for  ins- 
tance, 

Y*  = 0 and  Y**  = 1 » (65) 

Denote  by 

A*(t),  B*(t),  C*,  X*(t),o*,P*  (66) 

and 

A^^(t),  (t),  C**,  ^**(t)/  o**»  (67) 

the  particular  solutions  of  the  LTP-BVP  (48)- (56)  corresponding 
to  (65-1)  and  (65-2),  respectively.  Simple  manipulations , omitted 

for  the  sake  of  brevity,  show  that  the  general  solution  of 
(48)-(56),  valid  for  any  value  of  the  directional  coefficient 
Y,  can  be  written  as 

A(t)  = A* (t)  + Y [A** (t)  - A* (t) ] , (68-1) 

B(t)  = B*(t)  + YlB**(t)  - B*(t)  ] , (68-2) 

C = C*  + Y(C**  - C*)  , (68-3) 

and 

X(t)  = X*(t)  + YU**(t)  - X*(t)]  , (69-1) 
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+ 7(0**  - o*)  , 

(69-2) 

+ Y(y**  - u*)  . 

(69-3) 

As  a conclusion,  the  general  solution  of  (48) -(56)  requires 
that  two  sets  of  q + 1 sv;eeps  be  executed,  one  leading 
to  the  particular  solution  (66)  and  one  leading  to  the  parti- 
cular solution  (67) . 

5.7.  Stepsize  and  Directional  Coefficient.  With  the 
functions  A(t) , B(t),  C known,  the  following  two-parameter 
family  of  varied  functions  can  be  formed; 

x(t)  = x(t)+a|A*(t)  + Y[A**(t)  - A*(t)]  J , (70-1) 

u(t)  =u(t)+a|B*(t)  +Y[B**(t)  -B*(t)]|  , (70-2) 

IT  = ir  + a[C*  + Y(C**  - C*)  ] . (70-3) 

On  the  other  hand,  the  multipliers  X (t) , o , p form  the 
one-paraneter  family  (69)  . Upon  using  (69)  and  (70)  , we  see 
that  the  augmented  functional  (7)  takes  the  form 

J = J(a,Y)  . (71) 

Therefore,  the  optimum  values  of  a and  y satisfy  the  relations 


J^(a,Y)  =0, 


Jyio.,y)  =0  . 


(72) 


1 


u 


f 


i: 

[j 

0 

c 

[: 

[■ 

[■ 

c 

[ 

i 


27 


AAR-144 


Since  the  simultaneous  determination  of  a and  y might  be 
be  expensive  computationally,  we  proceed  in  a different  way. 
First,  we  determine  an  approximate  value  of  the  directional 
coefficient  y,  based  on  the  consideration  of  the  linear- 
quadratic  model  (Section  6).  Once  y is  known,  the  two-parameter 
feunily  (71)  reduces  to  the  one-parameter  family 

J*  J(a)  . (73) 

Then,  the  optimum  stepsize  a satisfies  the  relation 

J^(a)=0,  (74) 

whose  numerical  solution  can  be  obtained  using  quadratic  in- 
terpolation or  cubic  interpolation  (Ref.  9) . 


28 


AAR-144 


6.  Conjugate  Gradient  Phase;  Linear-Quadratic  Case 

In  the  previous  section,  we  analyzed  the  conjugate 
gradient  phase  in  general,  regardless  of  the  analytical  form 
of  the  functional  (1)  and  the  constraints  (2) -(5).  In  this 
section,  we  consider  the  linear-quadratic  case,  that  is,  the 
case  where  the  functional  (1)  is  quadratic  and  the  constraints 
(2) -(5)  are  linear. 

6.1.  General  Solution.  Under  the  assumption  of  linear 
constraints,  it  can  be  verified  that  the  particular  solutions 
(66)  and  (67)  satisfy  the  relations 

A**(t)-A*{t)  = A(t),  B**  (t)-B*  (t)=  B(t)  , C**-C*=C,  (75) 

^ t )—  X^(t)  = 0,  *^**~®*"®»  l^**~^*~®*  (76) 

As' a consequence,  Eqs.  (68)-(69)  reduce  to 

A(t)  = A*  (t)+ YA(t)  , B(t)=  B*  (t)+ YB(t)  , C = C*+ yC  , (77) 

X(t)=X*(t),  0 = 0*,  Vi  = Vi*.  (78) 

This  means  that  the  general  solution  of  the  LTP-BVP  (48) -(56) 
can  be  obtained  by  executing  only  one  set  of  q + 1 sweeps, 
namely,  the  set  of  sweeps  leading  to  the  solution  (66) . By 
the  way,  this  is  the  solution  corresponding  to  (65-1),  namely, 
the  solution  associated  with  the  ordinary  gradient  phase  of 


Ref.  5. 
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6.2.  Local  Orthogonality  Conditions.  Under  the 
assumption  of  linear  constraints,  the  two-parameter  family  (70) 
simplifies  to 


x(t)  =x(t)  +a[A*(t)  +YA(t)]  , 

(79-1) 

u(t)  =u(t)  +alB*(t)  +YB(t)]  , 

(79-2) 

if  = TV  + a (C*  + yC)  . 

(79-3) 

On  the  other  hand,  the  multipliers  X (t) , o,  y are  given  by 
Eqs.  (78).  Upon  using  (78)  and  (79),  we  see  that  the  augmented 
functional  (7)  still  takes  the  form  (71) . Hence,  the  optimum 
values  of  a and  y still  satisfy  the  relations  (72) . 

After  laborious  manipulations,  omitted  for  the  sake  of 
brevity,  Eqs.  (72)  lead  to  the  following  local  orthogonality 
conditions; 

1 

B*Bdt + cjc + (eJe)q = 0,  (80-1) 

0 

SB*Bdt  + cjc  + (eJe) 0 = 0 » (30-2) 

0 

with  the  implication  that 

1 

B*B*dt  + C*C*+  (eJe*)  Q =0  . (80-3) 
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Here,  the  adjective  local  is  employed  to  mean  that  Eqs.  (80) 
involve  vectors  B(t),  C,  E(0)  which  are  solutions  of  (48)- (56) 
computed  for  the  present  iteration  and  the  previous  iteration; 
they  also  involve  vectors  B* (t) , C*,  E*(0)  which  are  solutions 
of  (48)-(56)  for  y = 0 computed  for  the  present  iteration  and 
the  next  iteration. 

6.3.  Local  Conjugacy  Condition.  Let  w(t)  and  M(t)  denote 
the  vectors 


’x(t) 

A(t)  ' 

w(t)  = 

U (t) 

, M(t)  = 

B(t) 

TT 

C 

Let  f h denote  the  Hessian  matrices  of  the  functions 

ww  ww  ww 

f,g,h  with  respect  to  the  vector  w.  With  this  notation,  and 
under  the  assxamption  of  linear  constraints  and  quadratic  func- 
tional , Eqs.  (72)  lead  to  the  following  local  conjugacy  con- 
dition : 

SM'^f  Mdt  + (M'^h^Pl).  + (n'^g^M),  = 0 . (82) 

ww  ww  0 ww  i 

0 

Here,  the  adjective  local  is  employed  to  mean  that  Eq,  (82) 
involves  vectors  M(t) , that  is,  vectors  A(t) , B(t),  C, which 
are  solutions  of  (48) -(56)  computed  for  the  present  iteration 
and  the  previous  iteration. 
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6.4.  Stepsize  and  Directional ^Coefficient.  After 


observing  that 


M(t)  =M*(t)  + YM(t)  , 


and  after  accounting  for  Eqs.  (80)  and  (82)  , it  can  be  shown 
that  the  optimal  values  of  y and  a are  given  by 


where 


Y = Q/Q  , a = Q/R  , 


BjB*dt  + cJc*+(EjE*)Q  , 


(85-1) 


^ I AfP/N 

Q=  \ B*B*dt  + cic*+  (eX)o  ' 


(85-2) 


R = \ M^f  Mdt  + (m\  M )rt  + (M'^g  M),  , 

1 ww  ' ww  0 ^ww  1 


(85-3) 


Clearly,  the  optimal  directional  coefficient  y is  the 
ratio  of  the  error  in  the  optimality  conditions  Q for  the 
present  conjugate  gradient  iteration  to  the  error  in  the 
optimality  conditions  5 for  the  previous  conjugate  gradient 
iteration.  These  quantities  are  Icnown,  since  they  involve 
vectors  B*(t),  C*,  E*(0)  which  are  solutions  of  (48)- (56)  for 
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7 = 0 computed  for  the  present  iteration  and  the  previous  j 

iteration. 

, i 

With  7 known,  the  vector  M(t)  is  computed  with  (83)  and 
the  scalar  quantity  R is  computed  with  (85-3) . Then,  the 
optimal  stepsize  a is  evaluated  with  (84-2) . Clearly,  the 

optimal  stepsize  a is  the  ratio  of  the  error  in  the  optimality 

conditions  Q to  the  scalar  quantity  R,  which  constitutes  a 
measure  of  the  curvature  of  the  functional  (1) . Both  Q and  | 

R are  computed  employing  quantities  pertaining  to  the  present  ■ ■ 

conjugate  gradient  iteration. 

6.5.  Descent  Property.  Under  the  assumption  of  linear  | 

constraints,  Eq.  (60)  simplifies  to 

z = I B;Bdt  + c;c+  (e;e)q  . (86) 

Jo  . j 

Because  of  the  local  orthogonality  condition  (80-1)  written  | 

for  the  previous  iteration,  Eq.  (86)  yields 

Z=0.  (87) 


As  a consequence,  Eq.  (59)  reduces  to 

6J  = -oQ  , 


(88) 


where  the  error  in  the  optimality  conditions  Q is  given  by  Eq. 
(85-1).  Equation  (88)  holds  for  any  conjugate  gradient 
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iteration  and  shows  that,  since  Q>0,  we  have  6J  < 0.  Hence, 
for  a sufficiently  small,  it  is  guaranteed  that  the  augmented 
functional  J decreases.  In  conclusion,  for  the  linear- 
quadratic  case,  the  restart  procedure  mentioned  in  Section  5.4 
never  occurs.  This  means  that  the  directional  coefficient  y 
is  set  at  the  level  (61)  only  for  the  first  conjugate  gradient 
iteration. 

6.6.  General  Orthogonality  and  Conjugacy  Conditions. 

Now,  assume  that  the  algorithm  described  by  Eqs.  (48) -(56)  and 
(79)  is  employed,  starting  with  some  feasible  nominal  functions. 
Further,  assume  that  the  directional  coefficient  y is  set  at 
the  level  (61)  for  the  first  conjugate  gradient  iteration  and 
at  the  level  (84-1)  for  any  subsequent  conjugate  gradient 
iteration.  Under  these  assumptions  and  for  the  linear- 
quadratic  case,  one  can  generalize  the  local  orthogonality 
conditions  (80)  and  the  local  conjugacy  condition  (82)  as 
follows : 


0 

«1 


B*Bpdt  + C*Cp  + (E*Ep)  0=0/ 


^*^*p  0~  ® ' 


(89-1) 


(89-2) 


and 


M'^f  M dt+  )n  + (M'^g^M^),=  0,  (89-3) 

ww  p ww  p 0 ^ww  p 1 


0 
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where  the  subscript  p denotes  any  conjugate  gradient  iteration 
preceding  the  present  conjugate  gradient  iteration.  While 
these  equations  do  not  guarantee  convergence  in  a finite  num- 
ber of  steps,  they  do  guarantee  that  the  algorithm  generates 
a sequence  of  linearly  independent  vectors  M(t) , that  is,  a 
sequence  of  linearly  independent  variations  per  unit  stepsize 
A(t)  , B(t)  , C. 
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7 . Conjugate  Gradient  Phase;  Practical  Implementation 

In  this  section,  we  surnmarize  the  results  of  Sections 
5-6  and  suggest  practical  ways  of  utilizing  these  results, 
within  the  following  operating  rules:  (i)  the  use  of  second 
derivatives  must  be  avoided;  and  (ii)  the  simultaneous  de- 
termination of  the  directional  coefficient  y and  the  stepsize 
a must  be  avoided.  We  refer  to  the  general  case  where  the 
functional  (1)  is  nonquadratic  and/or  the  constraints  (2) -(5) 
are  nonlinear. 

7.1.  Auxiliary  Functions.  The  first  step  is  to  solve 
Eqs.  (48)- (56)  for  a fictitious  value  of  the  directional 
coefficient,  namely, 


Y*  = 0.  (90) 

Using  the  solution  technique  of  Section  4. 4, we  obtain  the 
following  auxiliary  functions  and  multipliers: 

A*(t),  B*(t),  C*,  A*(t),  o*  , y*  . (91) 

7.2.  Directional  Coefficient.  The  second  step  is  to 
compute  the  actual  value  of  the  directional  coefficient  y . 
For  the  first  conjugate  gradient  phase,  we  set 

Y = 0.  (92) 


For  subsequent  conjugate  gradient  phases,  we  set 
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Y = Q/Q  r 


(93) 


where 


Q = 


Q = 


BXdt  + cJc*-*-(EX)o  . 


ysHlA  aT*'^  A»PA 

+ C*C* + (EiE*)Q. 


(94-1) 


(94-2) 


In  Eqs.  (93) -(94),  the  symbols  Q and  Q denote  the  errors  in  the 
optimality  conditions  for  the  present  conjugate  gradient  phase 
and  the  previous  conjugate  gradient  phase,  respectively. 

The  directional  coefficient  (93)  is  acceptable  only  if 


J^(0)  = -(Q  + yZ)  < 0 , 

where  Q is  given  by  (94-1)  and  Z is  given  by 


(95) 


Z=\  B*Bdt  + cj6  + (eJe)  Q . (96) 

Jq 

If  Ineq.  (95)  is  violated,  then  the  directional  coefficient 
(93)  must  be  discarded  and  replaced  by  the  value  (92)  . This 
means  that  the  algorithm  must  be  restarted  by  replacing  the 
conjugate  gradient  phase  with  an  ordinary  gradient  phase. 


; ] 


[| 


7.3.  Basic  Functions.  The  third  step  is  to  compute  the 
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basic  functions  A(t),  B(t),  C and  the  nultipliers  X (t),  o,  u. 

This  is  done  with  the  following  formulas: 

A(t)  = A*(t)  + YA(t),  B(t)  = B*(t)  + YB(t),  C=C*  + yC,  (97) 

X(t)=X*(t),  a = a^,  ^ 

Therefore,  in  the  practical  implementation  of  the  algorithm, 
the  basic  functions  A(t),  B(t),  C are  computed  using  the  formulas 
derived  under  the  assumption  of  linear  constraints.  ; 

7.4.  Stepsize.  With  the  basic  functions  (97)  known,  we 
consider  the  one-parameter  family  of  varied  functions 

X ( t)  = X ( t)  + aA  (t) , u (t)  = u ( t)  + aB  (t),  tt  = u + aC  . (99) 

;l 

After  substitution  of  Eqs.  (98)-(99)  into  (7)  and  (14),  the 
following  functions  of  the  stepsize  are  obtained: 

J = J(a),  P = P(a).  (100) 

Then,  a one-dimensional  search  scheme  is  applied  to (100-1), 
and  a value  of  the  stepsize  a is  selected  for  which  the  fol- 
lowing relations  are  satisfied: 

J(«)<J(0),  P(a)<P*,  T(a)>0,  (101) 

where  t is  the  final  time  and  P*  is  a preselected  number,  not 
necessarily  small.  Satisfaction  of  Ineq.  (101-1)  is  possible 
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because  of  the  descent  property  of  the  conjugate  gradient 
phase.  Ineq.  (101-2)  is  introduced  to  prevent  excessive  cons- 
traint violation.  And  Ineq.  (101-3)  is  required  for  problems 
with  free  final  time. 

Prior  to  the  satisfaction  of  (101),  a scanning  process  is 
employed,  leading  to  the  bracketing  of  the  minimum  point  for 
J(a).  This  operation  is  then  followed  by  a Hermitian  cubic 
interpolation  process  (Ref.  9) , which  is  stopped  whenever  the 
following  relation  is  satisfied 


|j^(a)|<e3  or  I (a) /J„  (0)|  < , 


(102) 


subject  to  an  upper  limit  for  the  nvimber  of  search  steps  N . 

s 

Once  a stepsize  has  been  selected  consistently  with  either 
(102)  or  the  prescribed  upper  limit  for  the  number  of  search 
steps,  Ineqs.  (101)  must  be  checked.  If  satisfaction  occurs, 
then  the  stepsize  is  accepted.  If  any  violation  occurs, 

then  the  stepsize  must  be  bisected  progressively  until  sat- 

isfaction of  (101)  is  finally  achieved. 


The  symbols  and  denote  small,  preselected  numbers, 


3. 


Descent  Property  of  a Cycle 


A descent  property  exists  for  a complete  conjugate 
gradient-restoration  cycle  under  the  assumption  of  small 
stepsizes.  Let  denote  the  conjugate  gradient  stepsize  and 
the  restoration  stepsize.  Simple  manipulations,  omitted 
for  the  sake  of  brevity,  show  that  the  conjugate  gradient 
corrections  are  of  0(a  ),  while  the  restoration  corrections 

g 

2 

are  of  0{a^ag).  Hence,  for  sufficiently  small,  the  restor- 
ation corrections  are  negligible  with  respect  to  the  conjugate 
gradient  corrections.  Therefore,  the  restoration  phase  pre- 
serves the  descent  property  of  the  conjugate  gradient  phase. 

More  specifically,  let  I^,l2»  denote  the  values  of  the 
functional (1)  at  the  beginning  of  the  conjugate  gradient  phase, 
at  the  end  of  the  conjugate  gradient  phase,  and  at  the  end  of 
the  subsequent  restoration  phase.  Note  that  and  I2  are  not 
comparable,  since  the  constraints  are  not  satisfied  to  the 
same  accuracy.  On  the  other  hand,  and  are  comparable, 
and  the  conjugate  gradient  stepsize  can  be  selected  so 
that 


I3 < . (103) 

This  inequality  constitutes  the  descent  property  of  a complete 
conjugate  gradient-restoration  cycle.  In  order  to  enforce  it, 
one  proceeds  as  follows.  At  the  end  of  the  restoration  phase. 


1 
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one  must  verify  Ineq.  (103).  If  it  is  satisfied,  the  next 
conjugate  gradient  phase  is  started;  otherwise,  the  previous 


u 
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9.  Summary  of  the  Algorithm 

The  sequential  conjugate  gradient-restoration  algorithm 
solves  optimal  control  problems  involving  a functional,  sub- 
ject to  differential  constraints  and  general  boundary  condi- 
tions. The  algorithm  is  composed  of  a sequence  of  cycles, 
each  cycle  consisting  of  two  phases,  a conjugate  gradient 
phase  and  a restoration  phase.  The  objective  of  each  cycle 
is  to  decrease  the  functional  I,  while  the  constraints  are 
satisfied  to  the  predetermined  accuracy  (17-1) . 

The  decision  parameters  controlling  the  algorithm  are  the 
constraint  error  P and  the  optimality  condition  error  Q [see 
Eqs.  (14)  and  (15)].  If  P violates  Ineq.  (17-1),  the  algo- 
rithm executes  a restoration  phase.  If  P satisfies  Ineq. 

(17-1)  and  Q violates  Ineq.  (17-2) , the  algorithm  executes  a 
conjugate  gradient  phase.  Finally,  if  P and  Q satisfy  Ineqs. 
(17),  the  algorithm  stops:  convergence  has  been  achieved. 

9.1.  Restoration  Phase.  This  phase  involves  one  or 
more  iterations  and  can  be  summarized  as  follows. 

(a)  Assume  nominal  functions  x(t),  u(t),  ir  which  satis- 
fy condition  (3) , but  violate  at  least  one  of  conditions 

(2)  and  (4)-(5). 

(b)  For  the  nominal  functions,  solve  the  LTP-BVP 
(31)- (34)  and  (36) -(40)  using  the  method  of  particular  solu- 
tions. In  this  way,  obtain  the  functions  A(t),  B(t),  C and  the 
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multipliers  A(t) , o,  y. 

(c)  Using  the  functions  in  (b) , compute  the  restoration 
stepsize  by  a one-dimensional  search  on  the  constraint  error 
P(a).  To  this  effect,  perform  a bisection  process  on  a, 
starting  from  a^  = l,  until  Ineqs.  (47)  are  satisfied. 

(d)  Once  the  restoration  stepsize  is  known,  compute  the 
varied  functions  x(t),  u(t),  tt  with  Eqs.  (22). 

(e)  Verify  whether  the  varied  functions  in  (d)  satisfy 
Ineq.  (17-1) . If  this  is  the  case,  the  restoration  phase  is 
terminated.  Otherwise,  return  to  (a)  and  continue  the  pro- 
cess until  satisfaction  of  (17-1)  occurs. 

9.2.  Conjugate  Gradient  Phase.  This  phase  involves  a 
single  iteration  and  can  be  summarized  as  follows. 

(a)  Assume  nominal  functions  x(t),  u(t),  ir  which  satisfy 
the  constraints  (2) -(5)  within  the  preselected  accuracy  (17-1). 

(b)  For  the  nominal  functions  and  for  Y*  = 0,  solve  the 
LTP-BVP  (48) -(56)  using  the  method  of  particular  solutions. 

In  this  way,  obtain  the  auxiliary  functions  A*(t),  B*(t),  C* 
and  the  multipliers  X* (t) , a*,  y*  . 

(c)  Set  the  directional  coefficient  y at  the  level  (92) 
for  the  first  conjugate  gradient  phase  and  at  the  level  (93) 
for  any  subsequent  conjugate  gradient  phase . In  the  latter 
case,  accept  the  directional  coefficient  only  if  Ineq.  (95) 
is  satisfied;  otherwise,  reset  y at  the  level  (92). 
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(d)  Compute  the  basic  functions  A(t),  B(t),  C and  the 
multipliers  X (t) , a,  y using  Eqs.  (97)- (98). 

(e)  Using  the  functions  in  (d) , compute  the  conjugate 
gradient  stepsize  by  a one-dimensional  search  on  the  augmented 
functional  J (a ) until  satisfaction  of  Ineq.  (102)  occurs. 

Then,  bisect  the  resulting  stepsize  (if  necessary) , until 
satisfaction  of  Ineqs.  (101)  occurs. 

(f)  Once  the  conjugate  gradient  stepsize  is  known,  com- 
pute the  varied  functions  x(t),  u(t),  tt  with  Eqs.  (99). 

9.3.  Conjugate  Gradient- Restoration  Cycle.  After  the  restor- 
ation  phase  is  completed,  verify  whether  Ineq.  (103)  is  satisfied. 
If  this  is  the  case,  start  the  next  cycle  of  the  sequential 
conjugate  gradient-restoration  algorithm.  If  not,  return  to 
the  previous  conjugate  gradient  phase  and  reduce  the  conjugate 
gradient  stepsize  (using  a bisection  process)  until,  after 
restoration,  Ineq.  (103)  is  satisfied. 
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10 . Experimental  Conditions 

In  order  to  evaluate  the  theory,  nine  examples  were  solved . 
The  sequential  conjugate  gradient-restoration  algorithm  was 
programmed  in  FORTRAN  IV,  and  the  numerical  results  were  ob- 
tained in  double-precision  arithmetic. 

Computations  were  performed  at  Rice  University  using  an 
IBM  370/155  computer.  For  each  example,  the  interval  of  in- 
tegration was  divided  into  100  steps.  The  differential  equa- 
tions were  integrated  using  Hamming's  modified  predictor- 
corrector  method  with  a special  Runge-Kutta  starting  procedure 
(Ref.  10) . The  definite  integrals  I,  J,  P,  Q were  computed 
using  a modified  Simpson's  rule.  The  method  of  particular 
solutions  (Refs.  6-8)  was  used  to  solve  the  linear,  two-point 
boundary-value  problems  associated  with  both  the  conjugate 
gradient  phase  and  the  restoration  phase. 

10.1.  Convergence  Conditions . The  parameters  , £2^ 

12 

appearing  in  Ineqs.  (17)  and  (102)  were  set  at  the  levels 
ej^  = E-08,  e2  = E-04,  = E - 03  . (104) 

The  tolerance  level  (104-1)  characterizes  the  restoration 
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0 

0 

0 

I. 
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phase;  the  tolerance  levels  (104-1)  and  (104-2)  , employed  in 
combination,  characterize  the  algorithm  as  a whole;  and  the 
tolerance  level  (104-3)  characterizes  the  one-dimensional 
search  for  the  conjugate  gradient  stepsize. 

10.2.  Safeguards.  For  the  conjugate  gradient  phase, 
the  parameter  P*  appearing  in  Ineq.  (101-2)  was  set  at  the 
level 


P*  = 10.  (105) 

The  tolerance  level  (105)  limits  the  constraint  violation 
which  is  permissible  during  the  conjugate  gradient  phase. 

Also  for  the  conjugate  gradient  phase,  the  number  of  Hermitian 
search  steps  required  to  satisfy  Ineq.  (102)  was  subject  to 
the  upper  bound 


Ng  < 10  . (106) 

10.3.  Nonconverqence  Conditions.  The  sequential  conju- 
gate gradient-restoration  algorithm  was  programmed  to  stop 
whenever  satisfaction  of  any  of  the  following  inequalities 
occurred: 


(i) 

N > 50  , 

(107) 

(ii) 

> 30  , 

(108) 

1 


AAR-144 


(iii) 

Nj.  > 10  , 

(109) 

(iv) 

"br  " ' 

(110) 

(V) 

j^(0)  > 0, 

o 

II 

(111) 

(Vi) 

"bg  ’ “ ' 

Y=  0 , 

(112) 

(vii) 

“be  ^ ^0  ' 

Y = 0 . 

(113) 

Here,  N is  the  total  number  of  Iterations,  is  the  number  of 
cycles,  is  the  number  of  restorative  iterations  per  cycle, 
is  the  number  of  bisections  of  the  restoration  stepsize 
required  to  satisfy  Ineqs.  (47)  , is  the  number  of  bisec- 

tions of  the  conjugate  gradient  stepsize  required  to  satisfy 
Ineqs.  (101)  , is  the  number  of  bisections  of  the  conjugate 

gradient  stepsize  required  to  satisfy  Ineq.  (103),  and  Jqj(O) 
is  the  slope  of  the  augmented  functional  at  a=0. 

Inequalities  (107) -(108)  apply  to  the  algorithm  as  a 
whole.  Satisfaction  of  (107)  and/or  (108)  is  indicative  of 
extreme  slowness  of  convergence. 

Inequalities  (109) -(110)  apply  to  the  restoration  phase. 
Satisfaction  of  (109)  is  indicative  of  failure  to  produce  a 
feasible  solution  in  a reasonable  number  of  restorative  itera- 
tions. Satisfaction  of  (110)  is  indicative  of  extreme  small- 
ness of  the  restorative  displacements. 

Inequalities  (111) -(112)  apply  to  the  conjugate  gradient 


i ^ 1 


i • • 


47  AAR-144 

phase.  Satisfaction  of  (111)  means  that  the  descent  property 
of  the  conjugate  gradient  phase  does  not  hold,  owing  to  numeri- 
cal inaccuracy.  Satisfaction  of  (112)  is  indicative  of  extreme 
smallness  of  the  conjugate  gradient  displacements. 

Inequality  (113)  applies  to  a complete  conjugate 
gradient-restoration  cycle.  Satisfaction  of  (113)  is  indica- 
tive of  extreme  smallness  of  the  displacements  produced  within 
a complete  conjugate  gradient-restoration  cycle. 

10.4.  Restarting  Conditions.  The  directional  coefficient 
Y of  the  present  conjugate  gradient  phase  was  reset  at  the 
level  Y = 0 whenever  satisfaction  of  any  of  the  following  ine- 
qualities occurred; 


(i) 

J(,(0)  > 0 , 

Y = Q/Q  , 

(114) 

r~\ 

Y = Q/Q  , 

(115) 

LI  (iii) 

"be  ’ “ ' 

II 

o 

o> 

(116) 

0 Satisfaction  of 

(114)  means 

that  the  descent 

property  of 

1 r- 


the  conjugate  gradient  phase  does  not  hold.  Satisfaction  of 
(115)  is  indicative  of  extreme  smallness  of  the  conjugate 
gradient  displacements.  And  satisfaction  of  (116)  is  indica- 
tive of  extreme  smallness  of  the  displacements  produced  within 
a complete  conjugate  gradient-restoration  cycle. 


The  directional  coefficient  y of  the  next  conjugate 


n 
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..J 


gradient  phase  was  reset  at  the  level  y = 0 whenever  satisfac- 


tion  of  any 

of  the  following 

inequalities 

occurred : 

• 1 

(iv) 

1 < N.  < 10  , 

- bg- 

Y = 0 or 

y = Q/Q , 

(117) 

(V) 

- be  - ' 

Y * 0 or 

Y = Q/Q  . 

(118) 

■ f 

i : 


Satisfaction  of  (117)  or  (118)  is  indicative  of  large  viola- 
tions of  the  orthogonality  and  conjugacy  conditions,  owing  to 
the  fact  that  the  optimal  conjugate  gradient  stepsize  cannot 
be  employed. 
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11.  Nximerical  Examples 

In  this  section,  nine  numerical  examples  are  described 
employing  scalar  notation.  In  particular,  the  symbols  x^(t), 
i=l,...,n,  denote  the  components  of  the  state;  the  symbols 
u^(t),  i = l,...,m,  denote  the  components  of  the  control;  and 
the  symbols  i = l,...,p,  denote  the  components  of  the  para- 
meter. 

For  all  of  the  examples,  a time  normalization  is  used  in 
order  to  simplify  the  numerical  computations.  Specifically, 
the  actual  time  6 is  replaced  by  the  normalized  time 


t = e/T, 


(119) 


which  is  defined  in  such  a way  that  t = 0 at  the  initial  point 
and  t = 1 at  the  final  point.  The  actual  final  time  t,  if  it 
is  free,  is  regarded  as  a component  of  the  vector  parameter  tt 
to  be  optimized.  In  this  way,  an  optimal  control  problem  with 
variable  final  time  is  converted  into  an  optimal  control  prob- 
lem with  fixed  final  time. 

Example  11.1.  This  is  a linear-quadratic  problem  with 
(i)  initial  state  partially  given  and  (ii)  fixed  final  time 

T = 1: 


1 
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<1 


! 


( 


I 


1 


I 


\ ' ' 


i 

i 


Xi  = X2  + Ui  , X2  = Uj^  , 

(121) 

Xj^  (0)  = 2 , 

(122) 

X2  (1)  = 1 . 

(123) 

The  assumed  nominal  functions  are: 

Xj^(t)  = 2 , X2(t)  = 1 , Uj^(t)  = 

-1  . (124) 

The  numerical  results  are  given  In  Tables  1-2.  Convergence  to 

the  desired  stopping  condition  occurs  In  N = 4 Iterations, 

which  Include  1 restorative  Iteration  and  3 conjugate  gradient 

Iterations. 

Example  11.2.  This  Is  a linear-quadratic  problem  with 

(1)  a linear  relation  between  the  components  of  the  Initial 

state  and  (11)  fixed  final  time  t = 1: 

1=  (xj  + x^  + uj)dt  , 

(125) 

Xi  = X2  + Ui  , X2  = Uj^  , 

(126) 

Xj^  (0)  + X2  (0)  = 3 , 

(127) 

X2 (1)  * 1 . 

(128) 
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. li 


The  assumed  nominal  functions  are: 


(t)  = 2, 


X- (t)  = 1, 


Uj^  (t)  = -1. 


(129) 


The  numerical  results  are  given  in  Tables  3-4.  Convergence  to 
the  desired  stopping  condition  occurs  in  N = 4 iterations,  which 
include  1 restorative  iteration  and  3 conjugate  gradient  itera- 


tions. 


Example  11.3.  This  is  a problem  with  (i)  initial  state 
given  and  (ii)  fixed  final  time  t = 1: 


I = ^ (1  + x^  + X2  + u^)dt  , 


*1  " '^l  ■ *2  ' 


Xj^(O)  = 0, 


X2  = Ui  - x^X2  , 


X2(0)  = 1 , 


x^  (1)  =1,  X2  (1)  = 2 . 


The  assumed  nominal  functions  are: 


(130) 


(131) 


(132) 


(133) 


(t)  = t, 


X2  (t)  = 1 + t, 


Uj^(t)  = 1 . 


(134) 


The  numerical  results  are  given  in  Tables  5-6.  Convergence  to 
the  desired  stopping  condition  is  achieved  in  N = 7 iterations, 
which  include  5 restorative  iterations  and  2 conjugate 
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gradient  iterations. 

Example  11.4.  This  is  a problem  with  (i)  initial  state 
given  and  (ii)  fixed  final  time  t = 1; 

1=^  (-2cosu,)dt,  (135) 

Jo 


Xj^  = 2 sin  Uj^  - 1, 

*2  = , 

(136) 

X3^(0)  = 0, 

X 

to 

o 

II 

o 

(137) 

X 

I-* 

!-■ 

II 

o 

^2  (1)  =0.3  . 

(138) 

assumed  nominal  functions  are  : 

Xj^(t)  = 0,  *2^^^  “ 0.3t  , 

(t)  = 0 . 

(139) 

The  numerical  results  are  given  in  Tables  7-8.  Convergence  to 
the  desired  stopping  condition  occurs  in  N = 13  iterations, 
which  include  9 restorative  iterations  and  4 conjugate  gradient 
iterations. 

Example  11.5.  This  is  a problem  with  (i)  a nonlinear 
relation  between  the  components  of  the  initial  state  and  (ii) 
fixed  final  time  t = 1: 


= ^ (x^  ^2  ^ u^)dt  , 


*1  = *2  ^1  ' 


'2  ^1 


(140) 

(141) 


1! 


1.* 


L. 


0 
0 
0 

■ 

:[; 
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X2 (1)  = 1 . 


The  assumed  nominal  functions  are: 
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(142) 


(143) 


Xj^  (t)  = -2  , 


X2 (t)  = 1 , 


u^(t)  = 0 


(144) 


The  numerical  results  are  given  in  Tables  9-10.  Convergence 
to  the  desired  stopping  conditions  occurs  in  N = 7 iterations, 
which  include  4 restorative  iterations  and  3 conjugate  gradient 
iterations . 

Example  11.6.  This  is  a minimum  time  problem  with  (i)  a 
component  of  the  initial  state  given,  (ii)  a nonlinear  rela- 
tion between  the  remaining  components  of  the  initial  state, 
and  (iii)  free  final  time  t.  After  setting  = the  prob- 


lem is  as  follows: 


Xj^  — t 


I = , 


. 2 2 
X2  = (Xj^  - Uj^)  , 


X3  = TTj^  (Uj^-  X2  + Xj^), 


X^(0)  = 0, 


X2(0)  + X3(0)  = 1, 


X2(l)  = 0,  X3(1)=2, 


The  assumed  nominal  functions  are: 


(145) 


(146) 


(147) 


(148) 


'i-:K 
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Xj^(t)=0,  X2(t)=0,  X3tt)  = l+t,  UjCt)*!,  TTj^=l.  (149) 

The  numerical  results  are  given  in  Tables  11-12.  Convergence 
to  the  desired  stopping  condition  occurs  in  N=9  iterations, 
which  include  7 restorative  iterations  and  2 conjugate  gra- 
dient iterations. 

Excumple  11.7.  This  is  a minimum  time  problem  with  (i)  a 
component  of  the  initial  state  given,  (ii)  a nonlinear  relation 
between  the  remaining  components  of  the  initial  state,  and 
(iii)  free  final  time  t.  After  setting  the  problem  is 

as  follows: 


I = w^,  (150) 

* • • 

*1  ~ ^1^3°°®'^!'  X2  = TTj^x^sinUj^,  x^  = TTj^sinUj^  , (151) 

Xj^(0)=0,  X2(0)X3(0)=  0 , (152) 

x^(l)  = 1.  (153) 

The  assumed  nominal  functions  are: 

Xj^(t)=t,  X2(t)=l,  X3(t)=0,  U3(t)=l,  TTj^sl.  (154) 


The  numerical  results  are  given  in  Tables  13-14.  Convergence 
to  the  desired  stopping  condition  occurs  in  N=13  iterations, 
which  include  10  restorative  iterations  and  3 conjugate  gra- 
dient iterations. 


t 
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i. 

Example  11.8.  This  is  a minimum  time  problem  with  (i)  a 
component  of  the  initial  state  given,  (ii)  a nonlinear  rela- 
tion between  the  remaining  components  of  the  initial  state, 
and  (iii)  free  final  time  t.  After  setting  = the  prob- 


lem is  as  follows: 

I = TT^,  (155) 

Xl  = TilXsCosUi,  X2  = ^ j^x^sinuj^ , j^sinuj^ , (156) 

Xj^(O)  = 0,  X2(0)X3(0)  = 0,  (157) 

xjd)  + X3(l)  = 1 . (158) 


The  assumed  nominal  functions  are : 

Xj^(t)=t,  X2(t)=l,  X3(t)=0,  U3^(t)=l,  Tij^=l.  (159) 

The  numerical  results  are  given  in  Tables  15-16.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 9 iterations, 
which  include  7 restorative  iterations  and  2 conjugate  gradient 
iterations. 

Example  11.9.  This  is  a minimum  time  problem  with  (i)  a 
component  of  the  initial  state  given,  (ii)  a nonlinear  rela- 
tion between  the  remaining  components  of  the  initial  state, 
and  (iii)  free  final  time  t.  After  setting  1:3^=  t,  the 


n 
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problem  is  as  follows: 


I = TTj^  , 


(160) 


*1  ^ "^l^l' 


2 2 


X2  = TTj^  (Xj^  - , X,  = TT^  (U^  - + Xt  ) , (161) 


Xj^  (0)  = 0, 


x^(0)  +X3(0)  = 1 , 


X3(1)X2(1)  = 0, 


X3(l)  = 2. 


(162) 


(163) 


The  assumed  nominal  functions  are: 

X3(t)=t,  X2(t)=0,  X3(t)=l  + t,  U3(t)=l,  773=1.  (164) 

The  numerical  results  are  given  in  Tables  17-18.  Convergence 
to  the  desired  stopping  conditions  occurs  in  N = 9 iterations, 
which  include  7 restorative  iterations  and  2 conjugate  gra- 
dient iterations. 
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Table  1.  Convergence  history,  Exanqile  11.1. 


“c 

Nr 

N 

Y 

P 

Q I 

0 

0 

0 

0 

O.lOE+01 

1 

0 

1 

1 

0.22E-27 

0.92E+01  7.72083 

2 

1 

0 

2 

Y=0 

0.16E-28 

0.15E+00  5.89315 

3 

1 

0 

3 

YT^O 

0.40E-26 

0.84E-03  5.83868 

4 

1 

0 

4 

Y?^0 

0.77E-26 

0.24E-06  5.83848 

Table 

2. 

Converged  solution.  Example  11.1. 

t 

*1 

*2 

'^1 

0.0 

2.0000 

1.0198 

-2.0336 

0.1 

1.9137 

0.8410 

-1.5512 

0.2 

1.8580 

0.7081 

-1.1115 

0.3 

1.8333 

0.6175 

-0.7059 

0.4 

1.8408 

0.5661 

-0.3264 

0.5 

1.8819 

0.5516 

0.0337 

0.6 

1.9586 

0.5724 

0.3815 

0.7 

2.0736 

0.6277 

0.7230 

0.8 

2.2299 

0.7170 

1.0646 

0.9 

2.4313 

0.8408 

1.4124 

1.0 

2.6822 

1.0000 

1.7729 

T = 

1.00000 




.U  OJ  to  f-*  o 
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Table  3.  Convergence  history,  Excunple  11.2. 


:j 


N, 


**r 

N 

Y 

P 

Q 

I 

1 

0 

0 

0 

O.lOE+01 

J 

0 

1 

1 

0.13E-27 

0.97E+01 

5.99629 

1 

0 

2 

Y=0 

0.61E-29 

0.76E-01 

4.23041 

■ f 

1 

0 

3 

Y?^0 

0.40E-27 

0.19E-03 

4.20772 

1 

0 

4 

y^o 

0.14E-26 

0.30E-07 

4.20767 

• i 

.li 


Table  4. 

Converged 

solution.  Example 

11.2. 

t 

^1 

^2 

0.0 

1.2097 

1.7902 

-2.4220 

0.1 

1.1554 

1.5684 

-2.0209 

0.2 

1.1192 

1.3848 

-1.6564 

0.3 

1.1013 

1.2361 

-1.3225 

0.4 

1.1022 

1.1195 

-1.0137 

0.5 

1.1227 

1.0327 

-0.7250 

0.6 

1.1641 

0.9740 

-0.4516 

0.7 

1.2278 

0.9420 

-0.1891 

0.8 

1.3153 

0.9359 

0.0665 

0.9 

1.4290 

0.9552 

0.3198 

1.0 

1.5713 

1.0000 

0.5750 

i 

. i 


I 


T = 1.00000 


ij 

0 
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Table  5.  Convergence  history,  Exeunple  11.3. 


M 

g 

N 

Y 

P 

Q 

I 

0 

0 

0 

0 

0.72E+01 

1 

0 

4 

4 

0.32E-10 

0.97E+00 

33.67701 

2 

1 

1 

6 

Y=0 

0.84E-13 

0.50E-02 

33.46606 

3 

1 

0 

7 

Y?^0 

0.95E-08 

0.14E-04 

33.46465 

Table 

6 . Converged 

solution.  Example 

11.3. 

t 

*1 

*2 

0.0 

0.0000 

1.0000 

-8.3441 

0.1 

-0.7864 

0.2775 

-6.3703 

0.2 

-1.3015 

-0.2371 

-3.8640 

0.3 

-1.5841 

-0.5630 

-1.4827 

0.4 

-1.6737 

-0.7172 

0.4715 

0.5 

-1.6001 

-0.7107 

1.9964 

0.6 

-1.3776 

-0.5433 

3.2543 

0.7 

-1.0074 

-0.2050 

4.4921 

0.8 

-0.4872 

0.3184 

6.0506 

0.9 

0.1809 

1.0418 

8.4966 

1.0 

1.0000 

2.0000 

13.0501 

T = 1.00000 


o 1-1  N ro  Tj<  m 


000  O.lOE+01 

044  0.17E-08  0.67E+00  -1.11665 

127  Y=0  0.17E-11  0.34E-01  -1.16519 

119  yjiO  0.61E-09  0.31E-02  -1.16918 

1 1 11  Y?^0  0.73E-12  0.37E-03  -1.16953 

1 1 13  Y?^0  0.44E-15  0.47E-04  -1.16958 

I 


Table  8.  Converged  solution.  Example  11.4. 


t 

^1 

^2 

'^l 

0.0 

0.0000 

0.0000 

1.3386 

0.1 

0.0937 

0.0047 

1.3031 

0.2 

0.1855 

0.0186 

1.2615 

0.3 

0.2741 

0.0417 

1.1995 

0.4 

0.3572 

0.0733 

1.1095 

0 . 5 

0.4304 

0.1128 

0.9775 

0.6 

0.4839 

0.1587 

0.7553 

0.7 

0.4924 

0.2081 

0.3688 

0.8 

0.4146 

0.2543 

-0.1523 

0.9 

0.2383 

0.2877 

-0.6112 

1.0 

0.0000 

0.3000 

-0.8962 

T = 1.00000 
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Table  9.  Convergence  history.  Example  11.5. 


c 

N 

g 

N 

r 

N 

y 

P 

Q 

I 

0 

0 

0 

0 

O.lOE+01 

1 

0 

1 

1 

O.OOE-20 

0.47E+00 

3.33333 

2 

1 

1 

3 

Y=0 

O.lOE-16 

0.44E-02 

3.24276 

3 

1 

1 

5 

Y?^0 

0.35E-13 

0.28E-02 

3.24043 

4 

1 

1 

7 

Y?^0 

0.84E-12 

0.46E-05 

3.23912 

Table  10. 

Converged 

solution.  Example 

11.5. 

t 

^1 

^2 

'^l 

0.0 

-2.0314 

0.9345 

0.7621 

0.1 

-1.8703 

0.9987 

0.5292 

0.2 

-1.7250 

1.0418 

0.3391 

0.3 

-1.5935 

1.0677 

0.1851 

0.4 

-1.4739 

1.0798 

0.0617 

0.5 

-1.3647 

1.0809 

-0.0352 

0.6 

-1.2643 

1.0736 

-0.1089 

0.7 

-1.1712 

1.0598 

-0.1619 

0.8 

-1.0842 

1.0418 

-0.1958 

0.9 

-1.0015 

1.0213 

-0.2120 

1.0 

-0.9218 

1.0000 

-0.2114 

T = 1.00000 


Table  11.  Convergence  history.  Example  11.6 


■'ft  <tL 


- i 


N 

g 

N 

r 

N 

Y 

P 

Q 

I 

0 

0 

0 

0 

0.20E+01 

1 

0 

3 

3 

0.75E-11 

0.28E+00 

1.01089 

2 

1 

3 

7 

Y=0 

0.51E-14 

O.lOE-02 

0.86486 

3 

1 

1 

9 

Y?^0 

0.15E-10 

0.12E-04 

0.86430 

Table  12. 

Converged 

solution,  Example  11.6. 

t 

^2 

''3  "l 

0.0 

0.0000 

0.6465 

0.7628 

1.6143 

0.1 

0.1306 

0.4493 

0.8734 

1.4180 

0.2 

0.2465 

0.2970 

0.9937 

1.2700 

0.3 

0.3510 

0.1783 

1.1192 

1.1527 

0.4 

0.4463 

0.0871 

1.2475 

1.0539 

0.5 

0.5334 

0.0199 

1.3768 

0.9643 

0.6 

0.6130 

-0.0248 

1.5059 

0.8772 

0.7 

0.6850 

-0.0483 

1.6340 

0.7876 

0.8 

0.7490 

-0.0512 

1.7598 

0.6921 

0.9 

0.8044 

-0.0345 

1.8822 

0.5890 

1.0 

0.8506 

0.0000 

2.0000 

0.4782 

T = = 0 . 

86430 
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Table  13.  Convergence  history.  Example  11.7. 


N 

r 

N 

Y 

P 

Q 

I 

0 

0 

0 

0 

0.17E+01 

1 

0 

5 

5 

0.44E-16 

0.25E+00 

1.83370 

2 

1 

2 

8 

Y=0 

0.29E-09 

0.42E-01 

1.78266 

3 

1 

2 

11 

Y?^0 

0.29E-13 

0.16E-02 

1.77286 

4 

1 

1 

13 

y^o 

0.15E-11 

0.24E-05 

1.77245 

Table  14. 

Converged  solution , Example 

11.7. 

t 

*1 

*2 

^3 

"l 

0.0 

0.0000 

1.0000 

0.0000 

1.5693 

0.1 

0.0016 

1.0155 

0.1765 

1.4125 

0.2 

0.0129 

1.0607 

0.3486 

1.2556 

0.3 

0.0426 

1.1311 

0.5121 

1.0985 

0.4 

0.0974 

1.2198 

0.6630 

0.9416 

0.5 

0.1819 

1.3180 

0.7975 

0.7848 

0.6 

0.2974 

1.4163 

0.9125 

0.6281 

0.7 

0.4426 

1.5050 

1.0050 

0.4715 

0.8 

0.6129 

1.5755 

1.0728 

0.3147 

0.9 

0.8016 

1.6208 

1.1143 

0.1572 

1.0 

1.0000 

1.6363 

1.1280 

-0.0014 

T = = 1.77245 
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Table  15.  Convergence  history.  Example  11.8. 


N 

c 

N 

g 

Nr 

N 

Y 

p 

Q 

I 

0 

0 

0 

0 

0.17E+01 

1 

0 

4 

4 

0.15E-13 

0.45E-01 

1.02941 

2 

1 

2 

7 

Y=0 

0.25E-17 

0.32E-02 

1.00403 

3 

1 

1 

9 

y^o 

0.14E-08 

0.19E-04 

1.00002 

Table  16. 

Converged  solution,  Example  11.8. 

t 

*1 

*2 

*3 

"l 

0.0 

0.0000 

1.0000 

0.0000 

1.5648 

0.1 

0.0000 

1.0050 

0.1000 

1.5639 

0.2 

0.0001 

1.0200 

0.1999 

1.5630 

0.3 

0.0003 

1.0450 

0.2999 

1.5623 

0.4 

0.0006 

1.0800 

0.3999 

1.5615 

0.5 

0.0010 

1.1250 

0.4999 

1.5608 

0.6 

0.0016 

1.1800 

0.5999 

1.5600 

0.7 

0.0023 

1.2450 

0.6999 

1.5590 

0.8 

0.0033 

1.3200 

0.7999 

1.5579 

0.9 

0.0044 

1.4050 

0.8999 

1.5567 

1.0 

0.0058 

1.5000 

0.9999 

1.5551 

T = = 1.00002 

I 


I 


1. 


i I 


0 

f: 
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Table  17.  Convergence  history.  Example  11.9. 


Nc 

Nr 

N 

Y P 

Q 

I 

H 

H 

0 

0 

0.86E+00 

1 

0 

3 

3 

0.53E-10 

0.31E+00 

1.02610 

2 

1 

3 

7 

Y=0  0.24E-12 

0.69E-03 

0.86466 

3 

1 

1 

9 

Y?^0  0.33E-11 

0.48E-05 

0.86430 

Table  18. 

Converged  solution.  Example  11 

.9. 

t 

^1 

^2 

^3 

0.0 

0.0000 

0.6492 

0.7605 

1.6273 

0.1 

0.1314 

0.4496 

0.8718 

1.4247 

0.2 

0.2477 

0.2961 

0.9926 

1.2734 

0.3 

0.3524 

0.1770 

1.1185 

1.1545 

0.4 

0.4478 

0.0856 

1.2470 

1.0547 

0.5 

0.5351 

0.0185 

1.3766 

0.9645 

0.6 

0.6146 

-0.0260 

1.5059 

0.8768 

0.7 

0.6866 

-0.0492 

1.6339 

0.7868 

0.8 

0.7505 

-0.0518 

1.7598 

0.6910 

0.9 

0.8058 

-0.0349 

1.8822 

0.5876 

1.0 

0.8518 

0.0000 

2.0000 

0.4763 
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12.  Discussion  .j 

The  examples  presented  in  Section  11  were  solved  with  " ’ 

both  the  sequential  conjugate  gradient-restoration  algorithm 
(SCGRA)  of  this  report  and  the  sequential  ordinary  gradient- 
restoration  algorithm  (SOGRA)  of  Ref.  5.  This  was  done  in 
order  to  gain  perspective  on  the  relative  merit  of  SCGRA 
vis-a-vis  SOGRA. 

13 

The  comparative  results  are  presented  in  Tables  19-23, 
where  the  number  of  iterations  N required  to  achieve  different 
tolerance  levels  for  the  error  in  the  optimality  conditions  Q 
is  given  for  a fixed  tolerance  level  of  the  constraint  error 
P<E-08.  Also  shown  in  the  tables  are  the  values  obtained  for 
the  objective  functional  I. 

Cumulative  results  for  the  nine  examples  investigated  are  j 

given  in  Table  24.  Here,  the  total  number  of  iterations  for 

’ i 

convergence  EN  is  presented  as  a function  of  the  tolerance 
level  chosen  for  the  error  in  the  optimality  conditions  Q,  for 
a fixed  tolerance  level  in  the  constraint  error  P<E-08.  In 
this  comparative  study,  tolerance  levels  in  the  range  Q < E-02 
to  Q < E-06  were  chosen  for  the  error  in  the  optimality  condi- 
tions . 


13 


li 

f' 


In  Tables  19-23,  the  symbol  LQ  stands  for  a linear-quadratic 
problem,  and  the  symbol  NLQ  stands  for  a nonlinear  and/or 
nonquadratic  problem. 
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From  Tables  19-24,  the  following  observations  can  be 

made  . 

(i)  There  is  no  saving  in  number  of  iterations  for 

Q < E-02  and  Q<E-03.  The  saving  is  7%  for  Q<E-04,  14%  for 
Q<E-05,  and  18%  for  Q<E-06.  Clearly,  the  relative  advantage 
of  SCGRA  with  respect  to  SOGRA  increases  by  imposing  a tighter 
tolerance  level  on  the  error  in  the  optimality  conditions. 

(ii)  On  the  average,  the  relative  advantage  of  SCGRA  with 
respect  to  SOGRA  is  greater  for  problems  with  free  endpoints 
than  for  problems  with  fixed  endpoints. 

(iii)  On  the  average,  the  relative  advantage  of  SCGRA 
with  respect  to  SOGRA  is  greater  for  linear-quadratic  problems 
than  for  nonlinear  and/or  nonquadratic  problems. 

It  must  be  noted  that  the  experiments  performed  show  that 
the  computer  time  per  iteration  is  roughly  the  same  for  SCGRA 
and  SOGRA.  Therefore,  the  conclusions  pertaining  to  savings 
in  number  of  iterations  also  apply  to  savings  in  computer  time. 


Table  20. 

Results 

for  P < E-08 

and  Q < 

E-03. 

Example 

Type 

SCGRA 

SOGRA 

N 

I 

N 

I 

11.1 

LQ 

3 

5.83868 

4 

5.83854 

11.2 

LQ 

3 

4.20772 

4 

4.20768 

11.3 

NLQ 

7 

33.46465 

7 

33.46484 

11.4 

NLQ 

11 

-1.16953 

11 

-1.16950 

11.5 

NLQ 

7 

3.23912 

7 

3.23965 

]1.6 

NLQ 

9 

0.86430 

9 

0.86431 

:.1.7 

NLQ 

13 

1.77245 

10 

1.77262 

11.8 

NLQ 

9 

1.00002 

9 

1.00064 

11.9 

NLQ 

7 

0.86466 

7 

0.86466 

a 


us 

0 

D 


I 


I. 


( 


I I 


) 


Table  21.  Results  for  E-08  and  Q <:  E-04. 


Example 

Type 

SCGRA 

SOGRA 

N 

I 

N 

I 

11.1 

LQ 

4 

5,83848 

5 

5.83848 

11.2 

LQ 

4 

4,20767 

4 

4.20768 

11.3 

NLQ 

7 

33.46465 

7 

33.46484 

11.4 

NLQ 

13 

-1.16958 

13 

-1.16964 

11.5 

NLQ 

7 

3.23912 

11 

3.23921 

11.6 

NLQ 

9 

0.86430 

9 

0.86431 

11.7 

NLQ 

13 

1.77245 

12 

1.77245 

11.8 

NLQ 

9 

1.00002 

11 

1.00012 

11.9 

NLQ 

9 

0.86430 

9 

0.86430 

Table  22. 

Results  for  P < E-08 

and  Q < 

E-05. 

SCGRA 

SOGRA 

N I 

N 

I 

11.1 

LQ 

4 

5.83848 

6 

5.83848 

11.2 

LQ 

4 

4.20767 

5 

4.20767 

11.3 

NLQ 

8 

33.46464 

8 

33.46483 

11.4 

NLQ 

14 

-1.16962 

15 

-1.16965 

11.5 

NLQ 

7 

3.23912 

14 

3.23917 

11.6 

NLQ 

10 

0.86426 

10 

0.86429 

11.7 

NLQ 

13 

1.77245 

12 

1.77245 

11.8 

NLQ 

11 

1.00000 

14 

0.99998 

11.9 

NLQ 

9 

0.86430 

9 

0.86430 
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Table  23.  Results  for  P < E-08  and  Q<E-06. 


Exeunple  Type 


SCGRA 


SOGRA 


1 11.1 

LQ 

4 

5.83848 

6 

5.83848 

; 11.2 

LQ 

4 

4.20767 

6 

4.20767 

11.3 

NLQ 

8 

33.46464 

8 

33.46483 

■ 11.4 

NLQ 

15 

-1.16962 

17 

-1.16965 

11.5 

NLQ 

8 

3.23912 

17 

3.23917 

' 11.6 

NLQ 

10 

0.86426 

11 

0.86429 

11.7 

NLQ 

14 

1.77245 

13 

1.77245 

; 11.8 

NLQ 

12 

0.99997 

16 

0.99997 

i 11.9 

NLQ 

10 

0.86428 

10 

0.86429 

Table  24.  Cvunulative  number  of  iterations 
for  convergence^  P < E-08. 


SCGRA 

SOGRA 

EN 

EN 

56 

55 

69 

68 

75 

81 

80 

93 

85 

104 

13.  Conclusions 


In  this  report,  a new  member  of  the  family  of  sequential 
gradient-restoration  algorithms  for  the  solution  of  optimal 
control  problems  is  presented.  This  is  an  algorithm  of  the 
conjugate  gradient  type  and  solves  the  problem  represented  by 
Eqs.  (l)-(5);  Minimize  a functional  subject  to  differential 
constraints  and  general  boundary  conditions. 

The  algorithm  presented  here  differs  from  those  of  Refs. 
3-4,  in  that  it  is  not  required  that  the  state  vector  be 
given  at  the  initial  point.  Instead,  the  initial  conditions 
can  be  absolutely  general.  In  analogy  with  Refs.  3-4,  the 
present  algorithm  is  capable  of  handling  general  final  condi- 
tions; therefore , it  is  suitable  for  the  solution  of  optimal 
control  problems  with  general  boundary  conditions. 

The  importance  of  the  present  algorithm  lies  in  that  many 
optimal  control  problems  either  arise  naturally  in  the  present 
format  or  can  be  brought  to  such  a format  by  means  of  suitable 
transformations  (see  Ref.  2).  Therefore,  a great  variety  of 
optimal  control  problems  can  be  handled,  as  it  is  shown  by  the 
numerical  examples  presented. 

Nine  numerical  examples  are  presented  to  illustrate  the 
performance  of  the  algorithm.  The  numerical  results  show  the 
feasibility  as  well  as  the  convergence  characteristics  of  the 
algorithm.  A comparative  analysis  of  the  sequential  conjugate 


72 


AAR-144 


r 


( 


gradient-restoration  algorithm  and  the  sequential  ordinary 
gradient-restoration  algorithm  shows  that  the  relative  advantage 
of  SCGRA  with  respect  to  SOGRA  increases  by  imposing  a tighter 
tolerance  level  on  the  error  in  the  optimality  conditions. 

In  summary,  the  new  member  of  the  family  of  sequential 
gradient-restoration  algorithms  described  here  has  the  follow- 
ing properties;  (i)  it  retains  the  robustness,  reliability, 
and  convergence  characterisitics  of  the  algorithms  discussed  in 
Refs.  3-4;  (ii)  it  is  able  to  handle  all  of  the  optimal  control 
problems  treated  in  Refs.  3-4;  and  (iii)  it  has  the  additional 
capability  of  handling  optimal  control  problems  with  general 
boundary  conditions. 


I 


\ 
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